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Abstract 



In this paper we present results of using parallel supercom- 
puters to simulate beam dynamics in next-generation high 
intensity ion linacs. Our approach uses a three-dimensional 
space charge calculation with six types of boundary con- 
ditions. The simulations use a hybrid approach involving 
transfer maps to treat externally applied fields (including 
rf cavities) and parallel particle-in-cell techniques to treat 
the space-charge fields. The large-scale simulation results 
presented here represent a three order of magnitude im- 
provement in simulation capability, in terms of problem 
size and speed of execution, compared with typical two- 
dimensional serial simulations. Specific examples will be 
presented, including simulation of the spallation neutron 
source (SNS) linac and the Low Energy Demonstrator Ac- 
celerator (LEDA) beam halo experiment. 

1 INTRODUCTION 

The high intensity of future accelerator-driven systems 
places stringent requirements on the allowed beam loss, 
since very small fractional losses at high energy can pro- 
duce unacceptably high levels of radioactivity. Previous 
studies suggest that the low density, large amplitude halo 
of the beam is a major issue for these systems [|l], |J |]]. 
Large-scale simulations are an important tool for exploring 
the beam dynamics, predicting the beam halo, and facilitat- 
ing design decisions aimed at controlling particle loss and 
meeting operational requirements. 

The most widely used model for simulating intense 
beams in ion rf linacs is represented by the Poisson- 
Vlasov equations. These equations are often solved us- 
ing a particle-in-cell (PIC) approach. In this paper we will 
describe a parallel simulation capability that combines the 
PIC method with techniques from magnetic optics, and we 
will present results of using parallel supercomputers to sim- 
ulate beam dynamics in high intensity ion rf linacs. 

2 PHYSICAL MODEL AND 
NUMERICAL METHODS 

In the PIC approach a number of simulation particles, 
called macroparticles, are used to solve (indirectly) the evo- 
lution equations and model the charged particle dynamics. 
The motion of individual particles in the absence of radia- 
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tion can be described by Hamilton's equations, 

dq dH dp dH ^ 

dt dp 1 dt <9g' 

where H(q,p,t) denotes the Hamiltonian of the system, 
and where q and p denote canonical coordinates and mo- 
menta, respectively. In the language of mappings we would 
say that there is a (generally nonlinear) map, M, corre- 
sponding to the Hamiltonian H, which maps initial phase 
space variables, £\ into final variables, , and we write 

C f =MC. (2) 

The potential in the Hamiltonian includes contributions 
from both the external fields and the space-charge fields. In 
the Poisson-Vlasov approach, discreteness effects are ne- 
glected and the space charge is represented by a smoothly 
varying mean field. Typically, the Hamiltonian can be writ- 
ten as a sum of two parts, H = H ext + H sc , which corre- 
spond to the external and space-charge contributions. Such 
a situation is ideally suited to multi-map symplectic split- 
operator methods [Q] . A second-order-accurate algorithm 
for a single step is given by 

M(r) = M 1 (t/2) M 2 (t) Mi(t/2) , (3) 

where r denotes the step size, M \ is the map correspond- 
ing to H ext and M2 is the map corresponding to H sc . This 
approach can be easily generalized to higher order accuracy 
using Yoshida's scheme if desired [||]. 

The electrostatic scalar potential generated by the 
charged particles is obtained by solving Poisson's equation 



V 2 *(r) = -p(r)/e . 



(4) 



where p is the charge density. We have developed a 
Fourier-based transformation and an eigenfunction expan- 
sion method to handle six different boundary conditions: 
(1) open in all three dimensions; (2) open transversely and 
periodic longitudinally; (3,4) round conducting pipe trans- 
versely and open or periodic longitudinally; (5,6) rectangu- 
lar conducting pipe transversely and open or periodic lon- 
gitudinally. A discussion of the numerical algorithms for 
solving the Poisson's equation with these different bound- 
ary conditions can be found in . 

The charge density p on the grid is obtained by using 
a volume- weighted linear interpolation scheme^, Af- 
ter the potential and electric field is found on the grid, the 
same scheme is used to interpolate the field at the particle 
locations. During the course of the simulation each step in- 
volves the following: transport of a numerical distribution 
of particles through a half step based on M 1 , solving Pois- 
son's equation based on the particle positions and perform- 
ing a space-charge "kick" M2, and performing transport 
through the remaining half of the step based on Mi. 



3 APPLICATIONS 

We have applied the above 3D parallel PIC approach to an 
early design of the SNS linac and to the proposed LEDA 
beam halo experiment. Our simulation of the SNS linac 
starts at the beginning of the DTL. The code advances par- 
ticles through drift spaces, quadrupole fields and RF gaps. 
The dynamics inside the gaps is computed using external 
fields calculated from the electromagnetic code SUPER- 
FISH ^ . A schematic plot of the SNS linac configuration 
used in this study is shown in Figure i [0. It consists of 
three types of RF structures: a DTL, a CCDTL, and a CCL. 
There are a total of 425 RF segments in the linac. Figure 2 
shows the rms transverse size (x rmsj y r ms) and the maxi- 
mum transverse extent (x max , ymax) of the bunched beam 
in the linac with one set of errors. We see that the maxi- 
mum particle amplitude is well-below the aperture size of 
the linac. This margin is needed to operate the linac safely 
and to avoid beam loss at the high energy end. The jump 
in rms beam size between the DTL and CCDTL at 20 MeV 
is due to a change of focusing period from 8/3Atol2/3Aat 
805 MHz. 
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are two-fold: first, to study beam halo formation and test 
our physical understanding of the phenomena, and second, 
to evaluate our computational models and assess their pre- 
dictive capability through a comparison of simulation and 
experiment. Fig. 3 gives a schematic plot of the layout of 
the experiment It consists of 52 alternating-focusing 
quadrupole magnets with a focusing period of 41.96 cm. 
The gradients of the first four quadrupole magnets can be 
adjusted to create a mismatch that excites the breathing 
mode or the quadrupole mode. The transverse beam profile 
will be measured using a beam-profile scanner. Fig. 4 and 
Fig. 5 present simulation results of the transverse beam size 
for the breathing mode and the quadrupole mode, plotted 
at the center of the drift spaces between quadrupole mag- 
nets, as a function of distance. The plots include both the 
rms beam size and the maximum particle extent in the sim- 
ulation. The physical parameters for the simulation were 
1=100 mA, E=6.7 MeV, and f=350 MHz. The simula- 
tion was performed using 100 million macroparticles with 
a 128x128x256 (x-y-z) space-charge grid. 
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Figure 3: LEDA halo experiment layout 



Figure 1 : The SNS linac configuration 
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Figure 4: Transverse beam size as a function of distance 
for the breathing mode in the LEDA halo experiment 



Figure 2: Transverse beam size as a function of kinetic en- 
ergy in the SNS linac 

In the LEDA beam halo experiment, a mismatched 
high-intensity proton beam will be propagated through a 
periodic focusing transport system and measurements will 
be made of the beam profile. The goals of the experiment 



From Fig. 4, the two transverse components of the 
breathing mode are in phase, while the quadrupole mode 
in Fig. 5 has the two components out of phase. Evidently, 
it will be possible in the experiment to clearly excite ei- 
ther of the two modes. Furthermore, the debunching of 
the beam will not significanly alter the structure of the os- 
cillations. Fig. 6 shows the accumulated one-dimensional 
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100 million particle simulation of the LEDA halo experiment. 
Cumulative distributions of quadrupole mode oscillations. 



Figure 5: Transverse beam size as a function of distance 
for the quadrupole mode in the LEDA halo experiment 




Figure 6: Accumulated density profile along x for the 
breathing mode and the quadrupole mode 



density profiles (along x) for the breathing mode and the 
quadrupole mode just after the magnet #49. The breath- 
ing mode is more peaked and has a larger extent than the 
quadrupole mode. Measurements will be taken at this loca- 
tion and will be compared with our simulations. The data 
in Fig. 6 are well-resolved over a range of about 6 decades. 

An important piece of information from a design stand- 
point is amount of charge beyond a specified radius or spa- 
tial location as a function of distance along the accelerator. 
This is shown graphically in Fig. 7 which shows the com- 
plement of the horizontal and vertical cumulative density 
profiles, at every step, when the quadrupole mode is ex- 
cited. In other words, the contours describe the fraction of 
charge that would be intercepted by a scraper placed at that 
transverse position. 

The above LEDA simulations used 100 million 
macroparticles and a 3D Poisson solver, and required only 
2 hours to execute on 256 processors. In contrast, beam 
dynamics simulations performed on serial computers typi- 
cally use 10,000 to 100,000 macroparticles and a 2D Pois- 
son solver. Even if the above large-scale calculations could 
be performed on a PC, they would require on the order of a 
month to complete. In conclusion, while small-scale simu- 
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Figure 7: Horizontal and vertical cumulative density pro- 
files of a quadrupole mode mismatch in the LEDA halo 
experiment. 

lations on serial computers are extremely valuable for rapid 
design and predicting rms properties, large-scale simula- 
tions are needed for high-resolution studies aimed at mak- 
ing quantitative predictions of the beam halo. 
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